
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.0
 *  @date    Thu Feb 14 14:06:00 EST 2013
 *  @see     LICENSE (MIT style license file).
 *
 *  @see www2.cs.cas.cz/mweb/download/publi/Ples2006.pdf
 *  @see www.math.iit.edu/~fass/477577_Chapter_12.pdf
 *  @see Handbook of Linear Algrbra, Chapter 45
 */

// U N D E R   D E V E L O P M E N T  - to be merged with SVDecomp.scala

package scalation.linalgebra

import math.signum

import scalation.linalgebra.Householder.house
import scalation.linalgebra.MatrixD.{eye, outer}
import scalation.math.Basic.oneIf

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** This class is used to compute the Singlar Value Decomposition (SVD) of
 *  matrix 'a', i.e., decompose matrix 'a' into 'uu * ss * vv.t'.
 *  @param a  the m-by-n matrix to decompose (requires m >= n)
 */
class SVD (a: MatrixD)
      extends Error
{
    private val m = a.dim1         // number of rows
    private val n = a.dim2         // number of columns

//  if (n > m) flaw ("constructor", "SVD requires m >= n")

    private val vv = new MatrixD (m, m)
    private val uu = new MatrixD (n, n)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Use the Householder Bidiagonalization Algorithm to compute unitary matrices
     *  uu and vv such that uu.t * a * vv = b where matrix b is bidiagonal (nonzero
     *  elements on the main diagonal and the diagonal abobe it).
     *  @see Matrix Computations: Algorithm 5.4.2 Householder Bidiagonalization
     */
    def bidiagonalization (): MatrixD =
    {
        for (j <- 0 until n) {

            val (v, b) = house (a(j until m, j))         // Householder (vector, scalar)
            val ii  = eye (m-j)                          // identity matrix
            val bvv = outer (v, v) * b
            a(j until m, j until n) = (ii - bvv) * a(j until m, j until n)
            for (i <- j+1 until m) a(i, j) = v(i-j)

            if (j < n-2) {
                val (u, c) = house (a(j, j+1 until n))   // Householder (vector, scalar)
                val ii = eye (n-j-1)                     // identity matrix
                val cuu = outer (u, u) * c
                a(j until m, j+1 until n) = a(j until m, j+1 until n) * (ii - cuu)
                for (k <- j+2 until n) a(j, k) = u(k-j-1)
            } // if

        } // for

        for (i <- 0 until m; j <- 0 until n if j != i && j != i+1) a(i, j) = 0.   // clean
        println ("b = " + a)
        println ("vv = " + vv)
        println ("uu = " + uu)
        a                               // return b (a is changed in-place)
    } // bidiagonalization

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Use the Golub-Kahan Bidiagonalization Algorithm to compute unitary matrices
     *  u and v such that u.t * a * v = b where matrix b is bidiagonal (nonzero
     *  elements on the main diagonal and the diagonal abobe it).  Solve a * v = u * b
     *  by computing column vectors u_k and v_k and diagonals c and d.  Use transposes
     *  of u and v since row access is more efficient than column access.
     *  Caveat: assumes bidiagonals elements are non-negative and need to add
     *  re-orthogonalization steps.
     *  @see web.eecs.utk.edu/~dongarra/etemplates/node198.html
     */
    def bidiagonalization2 (): Tuple2 [VectorD, VectorD] =
    {
        var c  = new VectorD (n-1)                   // above main diagonal (beta)
        var d  = new VectorD (n)                     // main diagonal (alpha)
        var ut = new MatrixD (n, m)                  // transpose of u 
        var vt = new MatrixD (n, m); vt(0, 0) = 1.   // transpose of v
        val at = a.t                                 // transpose of matrix a

        for (k <- 0 until n) {
            ut(k)  = a * vt(k); if (k > 0) ut(k) -= ut(k-1) * c(k-1)
            d(k)   = ut(k).norm
            ut(k) /= d(k)

            val k1 = k+1
            if (k1 < n) {
                vt(k1)  = at * ut(k) - vt(k) * d(k)
                c(k)    = vt(k1).norm
                vt(k1) /= c(k)
            } // if
        } // for

        (c, d)
    } // bidiagonalization2

} // SVD class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** This object is used to test the SVD class.
 *  bidiagonalization answer = [ (21.8, -.613), (12.8, 2.24, 0.) ]
 */
object SVDTest extends App
{
     val a1 = new MatrixD ((4, 3), 1.,  2.,  3.,
                                   4.,  5.,  6.,
                                   7.,  8.,  9.,
                                  10., 11., 12.)
     var svd = new SVD (a1)
     println ("a1 = " + a1)
     println ("b1 = " + svd.bidiagonalization ())

//   val a2 = new MatrixD ((4, 3), 1.,  2.,  3.,
//                                 4.,  5.,  6.,
//                                 7.,  8.,  9.,
//                                10., 11., 12.)
//   svd = new SVD (a2)
//   println ("b2 = " + svd.bidiagonalization2 ())

} // SVDTest object

